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Abstract 


Pf  (o,  t)  =  p{3f  e  [o,  t]  :  g(x,  f((),  t )  <  0  }  (l) 


Using  the  total  probability  theorem,  we  propose  a  method  to 
calculate  the  failure  rate  of  a  linear  vibratory  system  with  random 
parameters  excited  by  stationary  Gaussian  processes.  The  response  of 
such  a  system  is  non-stationary  because  of  the  randomness  of  the 
input  parameters.  A  space-filling  design,  such  as  optimal  symmetric 
Latin  hypercube  sampling  or  maximin,  is  first  used  to  sample  the 
input  parameter  space.  For  each  design  point,  the  output  process  is 
stationary  and  Gaussian.  We  present  two  approaches  to  calculate  the 
corresponding  conditional  probability  of  failure.  A  Kriging 
metamodel  is  then  created  between  the  input  parameters  and  the 
output  conditional  probabilities  allowing  us  to  estimate  the 
conditional  probabilities  for  any  set  of  input  parameters.  The  total 
probability  theorem  is  finally  applied  to  calculate  the  time-dependent 
probability  of  failure  and  the  failure  rate  of  the  dynamic  system.  The 
proposed  method  is  demonstrated  using  a  vibratory  system.  Our 
approach  can  be  easily  extended  to  non-stationary  Gaussian  input 
processes. 

Introduction 

The  response  of  a  vibratory  system  with  random  parameters 
excited  by  stationary  Gaussian  processes  is  a  non-stationary  random 
process.  A  time-dependent  reliability  analysis  is  thus,  needed  to 
calculate  the  probability  that  the  system  will  perform  its  intended 
function  successfully  for  a  specified  time. 

Reliability  is  an  important  engineering  requirement  for 
consistently  delivering  acceptable  product  performance  through  time. 
As  time  progresses,  the  product  may  fail  due  to  time-dependent 
operating  conditions  and  material  properties,  component  degradation, 
etc.  The  reliability  degradation  with  time  may  increase  the  lifecycle 
cost  due  to  potential  warranty  costs,  repairs  and  loss  of  market  share. 
In  this  article,  we  use  time-dependent  reliability  concepts  associated 
with  the  first-passage  of  non-repairable  systems.  Among  its  many 
applications,  the  time-dependent  reliability  concept  can  be  used  to 
reduce  the  lifecycle  cost  [1]  or  to  set  a  schedule  for  preventive 
condition-based  maintenance  [2] . 

The  time-dependent  probability  of  failure,  or  cumulative 
probability  of  failure  [1,  3],  is  defined  as 


where  the  limit  state  g(X,Z(t),t)  =  0  depends  on  the  vector 
X  =  \Xx  X2  •••  Xm]  of  m  input  random  variables,  the 
vector  F(t)  =  [7] (/)  F2{t)  •••  Fn(t )]  of  n  input  random 

processes.  Failure  occurs  if  g(  )<  0  at  any  time  t  G  [o,  t]  where 
T  is  the  time  of  interest. 

The  time-dependent  probability  of  failure  of  Eq.  (1)  can  be 
calculated  exactly  as 

Pf  M  =  1  -  (l  -  ^})exp  j-  U(^j  (2) 

where  P®  is  the  instantaneous  probability  of  failure  at  the  initial  time 
and  A(t)  is  the  failure  rate.  Eq.  (2)  indicates  that  Pj.  M  can  be 

calculated  if  A(t)  is  known  and  vice  versa. 

In  the  commonly  used  up-crossing  rate  approach,  the  failure  rate 
is  approximated  by  the  up-crossing  rate 


1S„  P 
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under  the  assumptions  that  the  probability  of  having  two  or  more 
out-crossings  in  [m  +  Ar]  is  negligible,  and  the  out-crossings  in 

| ' t,  t  +  A^]  are  statistically  independent  of  the  previous  out- 
crossings  in  [(3,  r] .  Eq.  (2)  is  then  used  to  estimate  P j  (o,  7*) . 


Monte  Carlo  simulation  (MCS)  can  accurately  estimate  the 
probability  of  failure  of  Eq.  (1)  but  it  is  computationally  prohibitive 
for  low  failure  of  probability  problems.  To  address  the  computational 
issue  of  MCS,  analytical  methods  have  been  developed  based  on  the 
out-crossing  rate  approach  which  was  first  introduced  by  Rice  [4] 
followed  by  extensive  studies  [3,  5-7].  The  PHI2  method  [3]  uses  two 
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successive  time-invariant  analyses  based  on  FORM,  and  the  binomial 
cumulative  distribution  to  calculate  the  probability  of  the  joint  event 
in  Eq.  (3).  A  Monte-Carlo  based  set  theory  approach  has  been  also 
proposed  [8]  using  a  similar  approach  with  the  PHI2  method. 
Analytical  studies  such  as  in  [9,  10]  have  shown  that  the  PHI2 -based 
approach  lacks  sufficient  accuracy  for  vibratory  systems.  Other 
analytical  approaches  have  been  however,  proposed  to  estimate  the 
time-dependent  probability  of  failure  with  sufficient  accuracy  [11, 

12]. 

The  limited  accuracy  of  the  out-crossing  rate  approach  has  been 
improved  by  by  solving  an  integral  equation  involving  v+ {t)  and 

v++(mJ,  the  joint  up-crossing  rate  between  times  t  and  ^  [13]. 

The  up-crossing  rate  v+(r)  is  defined  in  Eq.  (3)  and  the  joint  up- 
crossing  rate  v++  ( t ,  )  is  defined  as 


u 


lim  - 
A/  — »  0 


°  [g  (• ,  t )  >  0  n  g  (• ,  t  +  At )  <  0  n  g  {• ,  f  j  >  0  n  g  [• ,  f  +  At  j  <  0  J 
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where  Y(^)  =  Y2{t)  *  *  •  is  the  vector  of  output 

(response)  random  processes, 

F(0=fo(0  F2{t)  •••  Fjt)]  is  the  vector  of  input  (force) 

random  processes,  and  X  =  \Xl  X2  *  *  *  Xm  ]  includes  m 
input  random  variables  (see  Figure  1).  The  mass  [m],  damping  [c] 
and  stiffness  [&]  matrices  depend  on  X.  Although  our  approach  can 
handle  all  n  random  processes  in  Y(^)  and  F (f) ,  we  will  consider 
the  case  with  only  one  input  Ft{t\  1  <i<n  and  one  output 
Yj  (t ),  1  <  j  <fl  for  simplicity.  Figure  1  provides  a  schematic 
with  the  input-output  notation. 


»Y(t) 


Figure  1 .  Schematic  of  linear  vibratory  system. 


(4) 

This  approach  has  been  adopted  in  [10]. 

Among  the  simulation-based  methods,  a  MCS  approach  was 
proposed  in  [14]  to  estimate  the  time-dependent  failure  rate  over  the 
product  lifecycle  and  its  efficiency  was  improved  using  an 
importance  sampling  method  considering  a  decorrelation  length  [15] 
in  order  to  reduce  the  high  dimensionality  of  the  problem.  Subset 
simulation  [16]  has  been  also  developed  as  an  efficient  simulation 
method  for  computing  small  failure  probabilities  for  general 
reliability  problems.  Its  efficiency  comes  from  introducing 
appropriate  intermediate  failure  sub-domains  to  express  the  low 
probability  of  failure  as  a  product  of  larger  conditional  failure 
probabilities  which  are  estimated  with  much  less  computational 
effort.  An  extreme  value  method  has  also  been  proposed  [17]  using 
the  distribution  of  the  extreme  value  of  the  response.  Recently,  a 
time-dependent  reliability  analysis  was  proposed  [18]  using  the  total 
probability  theorem  and  the  concept  of  composite  limit  state. 

In  this  paper,  we  present  a  time-dependent  reliability  analysis  for 
dynamic  systems  with  random  parameters  excited  by  a  stationary 
Gaussian  process  using  the  total  probability  theorem.  Metamodels  are 
used  to  estimate  conditional  probabilities  needed  by  the  total 
probability  theorem.  An  advantage  of  our  approach  is  that  we  can 
easily  handle  non-normal  and  correlated  random  variables  without 
additional  computational  effort. 

The  paper  is  arranged  as  follows.  Section  2  describes  the 
proposed  methodology  with  all  necessary  details.  Section  3  uses  a 
beam  example  to  demonstrate  all  developments  and  Section  4 
summarizes,  concludes  and  highlights  future  research. 


Proposed  Approach 

We  consider  the  following  n  degree-of-freedom  (DOF)  linear 
vibratory  system  with  random  parameters 


Total  Probability  Theorem  Approach 

The  calculation  of  the  time-dependent  probability  of  failure  (i.e., 
the  probability  of  the  response  exceeding  a  threshold)  is  very 
challenging  because  of  the  random  parameters  X.  For  each  realization 
of  X  however,  the  linear  vibratory  system  of  Eq.  (5)  has  constant 

coefficients  (matrices  [m],  [c]  and  [&]).  The  problem  then  becomes 
much  easier  if  the  system  is  excited  by  wide- sense  stationary 
processes  F (f)  .  The  total  probability  theorem  allows  us  to  calculate 
the  Pj  (o,r)  of  the  dynamic  system  with  random  parameters  in  terms 

of  conditional  probabilities  of  dynamic  systems  with  constant 
parameters. 

According  to  the  total  probability  theorem,  the  time-dependent 
probability  of  failure  of  Eq.  (1)  can  be  expressed  as  [18] 


Pf  (o,  t)  =  p(f)  =  \p(F/x).fx  (pdx 


(6) 


where  is  a  time-dependent  conditional  probability  of  failure, 

fx  (x)  is  the  joint  PDF  of  the  input  random  variables  X  and  Q  is 
the  support  of  /x(x) .  The  integral  of  Eq.  (6)  can  be  calculated 

using  numerical  integration  schemes  if  the  number  n  of  random 
variables  is  small  (e.g.  less  than  5).  Otherwise,  Monte  Carlo 
simulation  or  importance  sampling  methods  can  be  used.  In  all  cases, 
p(f/Q)  is  calculated  directly  or  using  a  pre-built  metamodel  as  is  the 

case  in  this  paper.  An  advantage  of  the  total  probability  theorem  is 
that  non-normal  and  correlated  random  variables  are  handled  without 
additional  computational  effort  or  loss  of  accuracy. 


[m(x)]Y(d+  [c(X)]Y(0+[qX)]Y(0  =  F  (f)  (5) 


Calculation  of  Conditional  Probability  p[f 


We  assume  that  the  input  random  process  F{t )  is  a  wide-sense 
stationary  and  Gaussian  random  process  with  zero  mean  and  a 
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constant  standard  deviation,  characterized  fully  by  the  power  spectral 
density  (PSD)  SFF  ( CD ) .  In  this  case,  for  each  realization  of  the 

input  random  variables  X,  the  output  process  Y{t )  is  also  wide- 

sense  stationary  and  Gaussian  with  zero  mean,  a  constant  standard 
deviation  and  a  PSD  of 

$yy  =  |  H  (^)|  $ FF  (ty)  (7) 

where  H{co )  is  the  frequency  response  function  of  the  system 
between  the  ith  (input)  and  jth  (output)  degrees  of  freedom  (DOF) 
[19]. 

For  a  real-life,  multi-degree  of  freedom  system  with  millions  of 
DOF,  the  calculation  of  H(co)  is  computationally  intensive.  It  is 

desired  therefore  to  minimize  the  number  of  times  H(cq)  is 


orthonormal  matrix  of  the  eigenvectors  O- ,/  =  1,. . N  and 
A  =  diag\Al  ^  •••  is  the  diagonal  matrix  of  the 

corresponding  eigenvalues.  Also,  let  Z  =  (zx  Z2  •  •  •  Z^)  be  a 

vector  of  N  independent  standard  normal  variables.  Because  of  the 
affine  transformation  property  of  the  multi-normal  distribution,  the 
following  spectral  representation  holds 

y(f)=//y(0+zVv<i>f(Tz,  (id 

i=i 

where  t  =  tl9 12,  . . tN  and  r<N  is  the  number  of  dominant 
eigenvalues.  Eq.  (6)  can  be  used  to  generate  sample  functions 
(trajectories)  of  Y(t).  The  Expansion  Optimal  Linear  Estimation 

method  (EOLE)  or  the  Orthogonal  Series  Expansion  (OSE)  [20]  can 
also  be  used. 


calculated.  Our  approach  using  the  total  probability  theorem  reduces 
the  number  H(cq)  is  calculated  considerably. 


According  to  the  Wiener-Khintchine  theorem,  the  Fourier 
transform  of  the  Syy{(jo)  spectrum  provides  the  autocorrelation 

function  /^(r)  of  the  output  process.  Because  both  5ry(ty)  and 
Ryy(t )  are  real  anc*  even  functions,  we  have  [19] 


+GO 

Ryy( t)  =  —  f  SYY{co)cos{coz)dT . 
n  i 


Note  that  for  a  wide-sense  stationary  process, 


(8) 


Since  Eq.  (11)  provides  a  discretized  version  of  the  output 
process,  we  can  use  Monte  Carlo  simulation  to  calculate  the  time- 
dependent  probability  of  failure.  However,  MCS  will  be 
computationally  very  demanding  considering  that  the  time  of  interest 
T  can  be  long  and  the  time  step  At  very  small.  To  reduce  the 
computational  effort,  we  use  in  this  paper  the  total  probability 
theorem. 


We  use  two  approaches  to  estimate  Ry/x)'  ^rst  °ne 

assumes  that  the  up-crossings  are  statistically  independent  while  the 
second  approach  does  not. 

Approach  1:  Statistically  Independent  Up-crossings 


RYY{z)=E\Y(t)-Y{t  +  Tt[  =  p<y+ /4  (9) 

where  —  1  < /2  < +1  is  the  correlation  coefficient 
and  £[r(r)]  =  JLly  and  <yY{t)  =  ° V  •  The  covariance  function  is 


CYY(t,t  +  —  RYY(t,t  +  /r)  —  +  r)]  — 

M  /  \  2  ^10) 

-  CYY  (t)  —  Ryy  (t)  —  /uY 

Knowing  the  covariance  function  of  the  output  process,  we  can 
use  a  spectral  decomposition  method  to  express  the  output  process  as 
a  linear  combination  of  the  eigenvectors  of  the  covariance  matrix 
where  the  coefficients  are  independent  standard  normal  random 

variables.  The  time  interval  of  interest  M  is  discretized  using  N 
discrete  times  tx  —  0,  £2?  •  •  •?  tN  =  T  and  the  NxN  covariance 
matrix  22  =  [c<9  v(/) ,  tj  ]J  ,  i  =  1,2, . . . ,  Af ;  j  =  1,2, . . Af  is 
formed  where  Cov(tt  ,tj  —  t(  +  T j  =  Cyy (tf ,  tj  j  =  Ryy (t )  is 
the  covariance  between  times  t  and  t  provided  by  Eq.  (10). 

1  J 


The  P\R/^J=  Pf[0,T]  is  estimated  using  Eq.  (2)  where  the 

failure  rate  (^)  is  approximated  by  the  up-crossing  rate  V+(/)  (i.e., 

x(t)  =  V+(/^) ).  Because  the  excitation  process  is  Gaussian  and  wide- 
sense  stationary  and  the  dynamic  system  is  linear  with  constant 
parameters,  the  up-crossing  rate  of  threshold  CC  is  constant  given  by 
[19] 


v+(t)  =  v+ 


where 


*  =  -2- 
2  7T(JV 


is  the  up-crossing  rate  for  CC  =  0  and 


<JY  =  J  J  co2SYY{oj)doj 


(12) 


(13) 


(14) 


Let  22  =  <I>  •  A  •  Or  be  the  spectral  decomposition  of  the 
covariance  matrix  22  where  0  =  [0X  02  •  •  •  0^  ]  is  the 


is  the  standard  deviation  of  the  derivative  process  Z(^) . 
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Approach  2:  Statistically  Dependent  Up-Crossings 

This  approach  demonstrates  the  significance  of  considering  the 
correlation  among  up-crossings  in  [0,  T\.  We  define  the  failure  event 


Because  the  integral  in  Eq.  (20)  is  non-negative,  fit)  has  the 
following  lower  bound 


as  F  =  {37  e  [o,  t]  :  g(x,  f(*  ),  t)  =  a -Y(t)<  0  } . 

The 

f(t)>v+(t)-f++(t,Tl)dr1.  (21) 

probability  [0,  t ],  0  <  t  <  T  is  then  given  by 

0 

■p1 

II 

+ 

1 

TsT 

(15) 

If  we  continue  the  process  of  using  Equations  (17)  and  (19)  to 
derive  Eq.  (20),  the  following  Rice’s  inclusion-exclusion  series  [4, 
19]  can  be  obtained 

where  f(t)  is  the  PDF  of  the  first  time  to  failure  and 

QQ  i~\  t  P\  1 

/(0  =  v+(0  +  ^](-l)  j  \--- \vl+(t,Tl,...,Ti_l)dTl---dzi_ 

'=2  00  0 

p‘J  =  p{f(o)>«} 

(16) 

(22) 

is  the  instantaneous  probability  of  failure  at  t  =  0 . 

The  up-crossing  rate  at  time  t  is  the  probability  that  the  first  up- 
crossing  occurred  at  t  or  at  a  previous  time  T  .  This  is  expressed  as 
[13,  10] 


where 


v'+Q  = 


v++0  'f  i  =  2 

vm0  if  i  =  3 


(23) 


v+(t)  =  f(t)  +  jv+(t\r)f(T)dr  (17) 

0 


The  PDF  fit)  is  then  used  in  Eq.  (15)  to  obtain  the  time-dependent 
probability  of  failure  P ^  (0,  t)  . 


where  v+  (t  I  t)  is  the  up-crossing  rate  at  t  conditioned  on  the  time 
of  the  first  up-crossing  being  T  .  Based  on  Eq.  (17),  f  (t)  has  the 
following  upper  bound 


f(t)<V+(t )  (18) 

because  the  integral  is  a  non-negative  quantity. 

The  joint  up-crossing  rate  V++(Y,Zj)  ,  indicating  the  probability 
that  an  up-crossing  occurs  at  both  t  and  Zj  <t ,  can  be  expressed  as 

Ti 

V++(t,T1)  =  V+(t\Tl)f(t)  +  ^V+(t,Tl\T2)f(T2)dT2  (19) 

0 

since  the  up-crossing  at  T1  is  the  first  or  the  first  up-crossing  has 

occurred  at  some  previous  time  T2  <  T\ .  Integration  of  Eq.  (19)  and 
substitution  in  Eq.  (17)  yields, 

t 

f(t)  =  V+(t)~\v++(t,Tl)dTl  + 

° 

t  T i 

+  |  ^V+(t,Tl\T2)f(T2)dT2dTl 

0  0 


For  the  series  of  Eq.  (22)  to  converge,  we  may  need  many  terms;  i.e., 
joint  up-crossing  rates  of  higher  order  (large  /).  To  avoid  this,  the 
integral  Eq.  (17)  can  be  modified  as 


t 


v+(t)  =  f(t)  +  j 


0 


v++(7,r) 

v+(r) 


f(T)dT 


(24) 


where  the  conditional  probability  definition  is  used  to  replace  the 


conditional  up-crossing  rate  v+  ( 1 1  t)  with 


v++(7,r) 

v+(r) 


where 


V+  (r)  is  the  up-crossing  rate  for  the  first  up-crossing.  Because  the 
up-crossing  rate  for  the  first  up-crossing  is  not  known,  an 
approximation  is  obtained  if  V+(r)  is  the  up-crossing  rate  for  any 


up-crossing  at  T  <t .  It  was  shown  in  [10,  13]  that  this 
approximation  provides  accurate  results.  The  example  in  this  paper 
verifies  this  claim. 


In  our  Approach  2,  we  solve  the  integral  Eq.  (24)  for  fit)  and 


then  use  it  in  Eq.  (15)  to  obtain  the  time-dependent  probability  of 
failure  P^  (0,  t)  . 


Note  that  if  all  up-crossings  are  assumed  statistically 
independent,  the  high  level  joint  up-crossing  rate  vl+  it ,  Tj , . .  -,Ti_l) 
in  Eq.  (22)  is  equal  to 
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vi+ (t,  Tx , . . r,_, )  =  V+  (/  )•  V+  (r,  )•  •  •  v+  (t,_,  ) .  (25a)  Vibratory  Beam  Example 


If  in  addition,  the  process  is  stationary,  the  up-crossing  rate  does  not 
depend  on  time  and 

v‘+(t,Tl,...,Ti_l)  =  v+(t}v+  (rt )•  •  •  v+  (r,_, )  =  (v+  )  .  (25b) 

In  this  case,  the  Rice’s  series  of  Eq.  (22)  yields 

fit)  =  v+  -  (v+)7  +  (v+]f -  (v+y  *—  + . . .  =  v+e~v+t .  (26) 

2  6 

Therefore,  the  failure  rate 


m=- 


nt ) 


v  e 


i-j  fwt 


\  +  [e~ 


-*°] 


=  V 


(27) 


is  constant  and  equal  to  the  up-crossing  rate  v+  . 

Calculation  of  v++(t,r)  for  Approach  2 


In  Eq.  (24),  we  need  the  up-crossing  rate  v+  (r)  and  the  joint 

up-crossing  rate  v++  ( t ,  T )  .  The  former  is  constant  for  each 
realization  of  X  and  is  calculated  using  Eq.  (12).  We  use  the 
following  steps  to  calculate  v++  ( t ,  r)  : 

•  Based  on  Eq.  (3),  P{St  O  Ft+At  )  =  V+At  where  St  indicates 
the  safe  event  at  t  and  Ft+At  indicates  the  failure  event  at 

t  +  At . 

•  Calculate  an  equivalent  reliability  index  ft  = —Q)  l(y+ At}. 
Because  of  stationarity  at  a  realization  of  X,  (3  is  the  same  for 
event  (SrnFT+At). 


Calculate  the  intersection  probability  Pt  T  using  the  bivariate 
standard  normal  CDF  02  as  [18] 

^=®2(-A-Ap)= 


1  ft  1  /  \ 

=  —  |  exp  — \J32  +  (32  -1(3  (3  cos^jcos ec20  d6 

arccos(/?)  -  ^ 

where  (3  is  obtained  from  the  previous  step  and  p  is  the 
correlation  coefficient  between  the  linear  safety  margins  at  t  and 
T  . 


We  demonstrate  the  proposed  approach  using  the  single  degree 
of  freedom  system  of  Figure  2  adopted  from  [21].  A  concentrated 
mass  M  is  placed  at  the  mid-span  of  a  massless  beam  of  length 
L  =  4m .  A  random  load  p(t )  is  applied  on  the  mass  which 

deflects  by  y(f)  as  shown  in  the  figure.  The  beam  is  made  of  steel 
with  Young  modulus  E  =  2\0GPa  and  density 


p  =  7850Kg/3  .  The  beam  has  a  rectangular  cross  section  of 
/  m 


width  h  and  height  h  and  provides  a  stiffness  K  = 


6EI 


to  the 


7  bti 

mass-damper  system  with  l  =  being  the  area  moment  of 

inertia.  There  is  also  a  damper  C  attached  to  the  mass. 


F(t) 


Figure  2  .  Beam  under  random  loading  [21] 

The  parameters  M,  b  and  h  are  normally  distributed  with 
pM  =100&g  ,  Pb  =  O.Olrn,  juh  =  0.05/77  and 
c JM  —  5  kg  ,  C7b  =  0.00  1  in ,  <yfi  —  0.005 m .  The  random  load 
F{t )  is  a  zero  mean,  wide-sense  stationary  process  with  the 
following  Pier son-Mo skowitz  spectrum 

$ff  (^)  =  e°J  (28) 

CO 

where  A  =  6.5  *109  N 2  sec-4  and  B  =  36,219sec  4  (see 
Figure  3). 


•  Based  on  the  definition  of  Eq.  (4),  we  have 

v++  0,r)  =  PlT(At)2 . 
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Figure  3.  Input  force  spectrum 

The  equation  of  motion  for  the  system  is 

My(t) + Cy{t)  +  Ky{t)  =  F(t) .  (29) 


The  undamped  natural  frequency  is  C0n  — 


A  damping  ratio 


g  =  0.3  is  used  resulting  in  a  damping  value  of  C  —  2 Mcong  .  The 

mean  value  of  the  undamped  natural  frequency  is  12.8  rad/sec.  The 
failure  event  is  defined  as 


F  =  {3f  e  [o,  T ] :  g(x,  F(/),  t)  =  0.02  -  Y(t )  <  0  } 


indicating  that  we  have  failure  if  the  response  y(f)  exceeds  the 
threshold  a  =  0.02  m . 


For  the  linear,  one  degree  of  freedom  vibratory  system  of  Eq. 
(29),  the  transfer  function  is 


H(cd) 


Y(<d) 


1 


F(co)  -Mco  +  jCco+K 

and  the  output  spectral  density  is 

S ff  {C0>) 


(30) 


SyyiCD)  =1  H(cd)  lz  Sff(co)  = 


(K-Mco2)2  +  CV 


(31) 


Figure  4  compares  the  conditional  time -dependent  probability  of 
failure  P^/^.j  f°r  M  =  100/:g  ,  b  —  0.01  m  and  h  =  0.05/71 

(mean  value  point  of  the  three  random  variables)  among  Approach  1, 
Approach  2  and  MC  simulation.  The  MC  results  have  been  obtained 
with  20,000,000  replications.  As  expected,  assuming  that  the  up- 
crossings  are  statistically  independent  and  using  Eq.  (2)  with 

49=v+  to  calculate  the  overestimates  the  actual 

probability.  However  the  results  of  Approach  2,  using  Eq.  (24)  and 
Eq.  (15),  are  very  close  to  MCS. 


Figure  4.  Comparison  of  conditional  time-dependent  probability  of  failure 
P^/^j  for  M-100  kg,  b=0.01  m  and  h=0.05  m 

Approach  1 


We  first  developed  a  Kriging  metamodel  for  V+  (Eq.  12)  in 
order  to  avoid  calculating  repeatedly  the  expensive  transfer  function 
H(co)  needed  to  calculate  the  output  spectrum  in  Eq.  (14). 

Note  that  although  H(co)  is  not  computationally  expensive  in  this 
example,  it  is  very  expensive  for  large-scale  vibratory  systems  with 
many  DOF.  We  also  developed  another  Kriging  metamodel  for  P^ 
in  Eq.  (2).  According  to  the  total  probability  theorem,  the  time- 
dependent  Pj  [o,r]  of  the  system  is  the  mean  value  of  the 

conditional  We  generated  10,000  replications  of  V+  and 

Pj  from  the  metamodels,  and  Eq.  (2)  was  used  for  each  replication 
to  calculate  p( y^-)  using  2(t)  ~  v+  . 

For  the  Kriging  metamodels,  we  used  an  Optimal  Symmetric 
Latin  Hypercube  (OSLH)  design  [22]  to  sample  the  design  space. 
The  range  of  values  for  each  random  variable  was 
\ju  —  3cr,  fd  +  3cr] .  We  started  with  20  design  points,  and  increased 
the  number  of  points  in  increments  of  10  until  convergence  of 
Pf\o,T  =  30]  was  reached.  The  “leave-one-out”  validation 

procedure  was  also  used  for  a  few  points  in  the  design  space  to  check 
the  accuracy  of  the  Kriging  metamodels.  Thirty  design  points  were 
needed  for  convergence  (see  Figure  5).  The  DACE  Matlab  Toolbox 
was  used  to  build  the  Kriging  metamodels.  A  second-order,  Gaussian 

correlation  structure  was  used  for  the  metamodel  of  V+ ,  and  a  first- 
order,  Gaussian  correlation  structure  was  used  for  the  metamodel  of 
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Figure  5.  Design  of  30  OSLH  points 

Figures  6  and  7  show  the  output  spectra  and  the  corresponding 
autocorrelation  functions  for  the  30  OSLH  design  points. 


Figure  6.  Output  spectra  for  the  30  OSLH  points 


Approach  2 

The  total  probability  theorem  of  Eq.  (6)  was  used.  To  determine 
the  time-dependent  conditional  probability  for  a  realization 

X  =  { M  b  h\  of  the  random  variables,  we  first  solved  the 
integral  Eq.  (24)  numerically  for  the  PDF  /(f)  of  the  first  time  to 
failure  for  0<f<30  using  a  Af  =  0.4sec  time  step. 
Subsequently,  Eq.  (15)  was  used  to  calculate 

p(/4)=//[°’4o<?<r. 

The  curve  was  calculated  for  each  point  of  an  OSLH 

design  similarly  to  Approach  1.  The  same  30  OSLH  points  were 
used.  Figure  8  shows  all  curves.  A  time-dependent  metamodel  [23] 
was  then  developed  using  a  singular  value  decomposition  method. 


Figure  8.  curves  for  the  30  OSLH  points  for  Approach  2 


Figure  9  compares  the  time-dependent  probability  of  failure 
between  the  two  approaches.  Because  of  the  statistically  independent 
assumption  of  the  up-crossings,  the  Approach  1  overestimates  the 
probability  of  failure.  The  results  from  Approach  2  are  much  more 
accurate  as  shown  in  Figure  4  where  Approach  2  is  much  closer  to 
MCS  for  the  case  where  all  random  variables  are  at  their  means. 
Figure  10  compares  the  time-dependent  failure  rate  calculated  using 
the  two  approaches.  Figure  11  shows  the  joint  up-crossing  rate 

V++(f,r)  for  the  design  point  with  all  random  variables  at  their 
mean  values. 


Figure  7.  Autocorrelation  function  for  the  30  OSLH  points 


Page  7  of  9 


DISTRIBUTION  STATEMENT  A.  Approved  for  public  release;  distribution  is  unlimited. 


1 


DISTRIBUTION  STATEMENT  A.  Approved  for  public  release;  distribution  is  unlimited. 


Figure  9.  Comparison  of  time-dependent  probability  of  failure 
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Figure  10.  Comparison  of  time-dependent  failure  rate 
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Figure  1 1 .  Joint  up-crossing  rate  for  design  point  with  all  random  variables  at 
their  means 


Summary,  Conclusions  and  Future  Work 

We  presented  a  methodology  to  calculate  the  time-dependent 
probability  of  failure  and  the  failure  rate  of  a  linear  vibratory  system 
with  random  parameters  excited  by  stationary  Gaussian  processes 
using  the  total  probability  theorem  and  an  integral  equation  involving 
the  up-crossing  and  joint  up-crossing  rates.  An  optimal  symmetric 
Latin  hypercube  space-filling  design  is  first  used  to  sample  the 
system  random  parameters.  Time-dependent  conditional  probabilities 
are  then  calculated  at  each  design  point  by  solving  an  integral 
equation  involving  up-crossing  and  joint  up-crossing  rates  as  well  as 
the  PDF  of  the  first  time  to  failure.  A  time-dependent  metamodel  of 
the  conditional  probabilities  is  built  and  used  in  the  total  probability 
theorem  to  calculate  efficiently  and  accurately  the  time-dependent 
probability  of  failure  and  the  failure  rate  of  the  vibratory  system. 

Our  approach  assumes  a  wide- sense  stationary  and  Gaussian 
excitation.  However,  it  can  be  easily  extended  to  handle  non- 
stationary  Gaussian  excitations.  An  example  of  a  simple  vibratory 
system  was  used  to  demonstrate  all  developments. 
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